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ABSTRACT 


12 A new dynamical core of Environment and Climate Change Canada’s 

13 Global Environmental Multiscale (GEM) atmospheric model is presented. 

14 Unlike the existing log-hydrostatic-pressure-type terrain-following vertical 
is coordinate, the proposed core adopts a height-based approach. The move to 

16 a height-based vertical coordinate is motivated by its potential for improving 

17 model stability over steep terrain, which is expected to become more prevalent 
is with the increasing demand for very high resolution forecasting systems. A 

19 dynamical core with height-based vertical coordinate generally requires an it- 

20 erative solution approach. In addition to a three-dimensional iterative solver, a 

21 simplified approach has been devised allowing the use of a direct solver for the 

22 new dynamical core that separates a three-dimensional elliptic boundary value 

23 problem into a set of two-dimensional independent Helmholtz problems. The 

24 new dynamical core is evaluated using numerical experiments that include 

25 two-dimensional nonhydro static theoretical cases as well as 25 -km resolution 

26 global forecasts. For a wide range of horizontal grid resolutions—from a 

27 few meters to up to 25 km—the results from the direct solution approach is 

28 found to be equivalent to the iterative approach for the new dynamical core. 

29 Furthermore, results from the numerical experiments confirm that the new 
so height-based dynamical core leads to results that are equivalent to the existing 
31 pressure-based core. 
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32 1 . Introduction 

33 The dynamical core of the Global Environmental Multiscale (GEM) model, used operationally 

34 by Environment and Climate Change Canada (ECCC) for numerical weather prediction (NWP), 

35 employs a log-hydrostatic-pressure-type terrain-following vertical coordinate. The system of 

36 nonlinear model equations is linearized around a basic state and is then reduced to an elliptic 

37 boundary value (EBV) problem through numerical discretization and elimination of variables 

38 (Girard et al. 2014 ). The existing pressure-type coordinate system then permits the use of a direct 

39 solver for the discretized EBV problem to resolve the dynamical component of the flow. The 

40 direct solver starts by separating the EBV problem vertically in terms of the vertical eigenvectors 

41 of the part of the coefficient matrix that only includes the discretized difference and average 

42 operators in the vertical direction (Qaddouri and Lee 2010 ). For N number of model vertical 

43 levels, the resulting N vertically-decoupled two-dimensional Helmholtz problems are then 

44 separated along the longitude, leading to a system of tridiagonal problems depending only on the 

45 latitude. The tridiagonal problems are finally solved using LU decomposition without pivoting. 

46 Such an approach is computationally more efficient than most iterative methods, particularly for 

47 the spatial resolutions of the current operational NWP systems at ECCC. 

48 

49 One of the principal incentives for the adoption of the existing pressure-type vertical coordinate 
so in GEM is the computational advantage of the direct solution approach that is permissible with 
si such a coordinate. However, numerical experiments carried out at ECCC have revealed that 

52 vertical separability, which is an imperative for the direct solution approach, can become quite 

53 restrictive for very high spatial resolutions, e.g., for sub-kilometer horizontal grid spacing. Fur- 

54 thermore, the reduction of the 2 D Helmholtz problems into the ID tridiagonal problems requires 
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ss projection of right hand side (RHS) of the 2 D problems along the pertinent eigenvectors which is 

56 based on Fourier transformation. This necessitates transposing the coefficient matrix that involves 

57 global communications, and therefore, becomes inefficient for very large number of processor 
ss cores. As a result, the direct solver is found to lose scalability with increasing number of processor 

59 cores. Initial research at ECCC as well as published research literature (Muller and Scheichl 2014 ) 

60 reveal that optimized three-dimensional iterative solvers may possess better scalability in these 

61 circumstances. The different limitations of the existing direct solver, particularly its potential lack 

62 of scalability for future generations of massively parallel supercomputers, therefore warrants the 

63 development of more scalable iterative solvers at ECCC. A height-based vertical coordinate is 

64 considered to be more amenable to such iterative solvers as the metric terms originating from the 
es vertical coordinate transformation appear explicitly in the discretized EBV. 

66 

e? Another, but more challenging, problem pertaining to the existing GEM dynamical core is 
es the fact that the current model exhibits strong numerical instability over steep orography. Tests 

69 conducted at ECCC have determined the maximum permissible terrain slope for maintaining 

70 stability to be approximately 45 ° (Vionnet et al. 2015 ). This is generally considered to be a 

71 limitation inherent to the terrain-following coordinate (TFC) systems (Zangl 2012 ). With a 

72 growing demand for very high-resolution operational NWP systems, steep orographic slopes 

73 are expected to become more prevalent in the near future. Improving model stability over steep 

74 mountains is therefore of critical importance for the future development of sub-kilometer NWP 

75 systems. A number of approaches have been investigated to improve numerical stability with the 

76 existing pressure-type vertical coordinate in GEM. These include increased off-centering in the 

77 discretized vertical momentum equation, a vertically-variable basic state temperature profile, and 

78 modifications to the nonhydro static contributions in the linear system arising from the discretized 
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79 GEM formulation. However, none of these approaches has been found to lead to any meaningful 
so stability improvement for steep orography. 

81 

so Although the instability induced by steep orography is often characterized as a limitation of 
as the terrain-following nature of the vertical coordinate itself, Smolarkiewicz et al. ( 2007 ) were 

84 successful in resolving flow around the Pentagon with a model involving height-based TFC 

85 where the maximum slope was well above the widely acknowledged 45 ° threshold. Previous 

86 experience with the Mesoscale Compressible Community model at ECCC also suggests that a 
a? dynamical core with height-based TFC does not suffer from similar severe orography-induced 
as instability (Girard et al. 2005 ). Apparently, a dynamical core with a height-based TFC can lead 

89 to improved numerical stability through better implicit treatment of the metric terms arising from 

90 the vertical coordinate transformation through iterative solvers. More importantly, conventional 

91 numerical approximation of the horizontal gradients in a TFC becomes less accurate with 

92 increasing terrain slope as well as with increasing vertical resolution close to the model surface 

93 (Mahrer 1984 ). In this context, Zangl ( 2012 ) argues that the pressure gradient term, in particular, 

94 becomes susceptible to triggering numerical instability when the height difference between two 

95 adjacent grid points along a terrain-following vertical level is much larger than the vertical grid 

96 resolution adjacent to the level. Numerical approximation of the horizontal gradient terms in 

97 the TFC, however, can be significantly improved following the corrections proposed by Mahrer 

98 ( 1984 ). These corrections require determination of the modified horizontal differencing stencils 

99 associated with each grid-point location that minimize the error in the metric corrections for the 

100 terrain-following nature of the coordinate. The existing pressure-type TFC varies with time and, 

101 therefore, would require determination of the pertinent grid-point locations in the vertical for 

102 the modified differencing stencils at every time step. On the contrary, the height-based TFC is 
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103 


time-invariant and thus would require the determination of these grid-point locations only once 

104 at the beginning of the time integration. Therefore, from a computational efficiency standpoint, 

105 a height-based TFC is more suitable for implementing improved numerical approximation of 

106 horizontal gradients to address instabilities induced by steep orography. 

107 

108 The aforementioned challenges associated with the existing log-hydrostatic-pressure-type 

109 TFC motivated the development of a new dynamical core for the GEM model that utilizes a 
no height-based TFC. The primary objective of the present study is to demonstrate that, for the 
m model configurations where orography-induced numerical instability is not relevant—i.e., for 
112 horizontal grid resolutions within the hydrostatic regime—the new dynamical core developed at 
ns ECCC with height-based TFC makes predictions that are equivalent to those from the existing 
n4 model. The present study also explores the appropriate strategy for coupling the operational 
ns Physics Parameterization Package (PPP) of RPN (Recherche en prevision numerique) with the 
ns new height-based dynamical core. Different setups for numerical experiments are utilized to 
n7 compare the newly-developed dynamical core with the existing one covering both hydrostatic and 
ns nonhydrostatic scenarios. The experiments include two-dimensional theoretical cases (Robert 
n9 1993 ; Schar et al. 2002 ) as well as three-dimensional global forecasts over a Yin-Yang grid 

120 (Qaddouri and Lee 2011 ). 

121 

122 Relevant background information on the GEM dynamical core with the proposed height-based 

123 TFC—from the spatial and temporal discretizations to the derivation of the elliptic boundary value 

124 problem—is presented in section 2 . The different solution methods utilized for the discretized 

1 25 elliptic problem is discussed in section 3 . The issue of coupling between the dynamical core and 

126 the parameterized physics forcings is presented in section 4 . Section 5 contains the comparisons 
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127 between the existing and the proposed dynamical cores in the context of two-dimensional 

128 theoretical benchmark cases as well as three-dimensional deterministic global predictions. The 

129 conclusions are then summarized in section 6. 


131 2 . Model Description 

132 a. Governing equations 

133 The GEM model equations originate from the Euler equations. With the traditional shallow at- 

134 mosphere approximation (Phillips 1966 ), the system of equations in a spherical coordinate (A, (j ), r) 

135 can be expressed as follows: 


136 
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( 1 ) 

( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 


mo where Eqs. (l)-( 5 ) govern the evolutions of the u, v, and w components of velocity, mass and 

141 energy, respectively. The spatial coordinates in the above equations are denoted by (x.y,z) which 

142 are related to the spherical coordinate (A, 0, r) through the differential relations given by 


dx — acos(j)d?i,dy = ad<f),dz = dr, 


( 6 ) 
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143 such that u, v and w are the physical wind components. In Eq. (6), a denotes Earth’s radius. The 

144 Lagrangian derivative in this case can be expressed as 


d d d d d 

— — ^ -b -bV^r-bW^—. 

dt dt dx dy dz 


( 7 ) 


In addition to the four independent variables (jc,y, z, t), the system of five equations (l)-( 5 ) involves 
six dependent variables, namely, the velocity components u, v, and w, the temperature T. the 
pressure p, and the density p. Also, in the above equations, / is the Coriolis parameter and 
K = R/cp where R is the gas constant and c p is the specific heat of air at constant pressure. The 
terms on the RHS of Eqs. (l)-( 5 ) with subscript “ phys ” denote the various physical forcings. 
Depending on the equation, these physical forcings may arise from different sources that include 
friction, diabatic heating and frictional dissipation of kinetic energy. A sixth equation is required 
to close the system described by the six dependent variables and is provided by the equation of 
state, given by 


p = pRT. 


( 8 ) 


154 It is important to note that the atmospheric substance does not only contain dry air but also water 

155 vapor and different types of hydrometeors. The displacement and evolution of water vapor and hy- 

1 56 drometeors in the atmosphere are governed by their own evolution equations. However, they will 
is? also affect the RHS terms through fluxes of water vapor and precipitation which constitute sources 
i 58 of mass. The total air density in the presence of water in its different forms can be expressed as 


P — Pd ~b Pw — Pd ( 1 ~b r w ), ( 9 ) 

159 where p c / is the dry air density, p w — ^ is the density of water vapor and hydrometeors, and r w 

160 is the mixing ratio for the total water content of the atmosphere. The equation of state in such a 
lei scenario is strictly given by 

P= (PdRd + p v R v )T, ( 10 ) 
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162 where p v and R v are the density and gas constant of water vapor. Eq. ( 10 ) can then be further 

1 63 rearranged in terms of the total air density p as 

p = pR d T v , ( 11 ) 

1 64 where T v is the virtual temperature of moist air which is given by 

1 + ^: r v 

T v = 1 Kd T. ( 12 ) 

1 65 where r v — ^ is the water vapor mixing ratio. Rewriting the dynamical equations (l)-( 5 ) in terms 
lee of T v is helpful as the equations can then be expressed using only the dry gas constant that does 

1 67 not vary due to the atmospheric water content. It also allows to account for the effects of water 

1 68 vapor buoyancy and condensed water loading implicitly. Furthermore, from the adiabatic point 

1 69 of view, the introduction of T v has no consequence since the water content is conserved during 

170 dynamical transport. 

171 

172 b. Vertical coordinate 

173 The log-hydrostatic-pressure-type terrain-following vertical coordinate of the operational GEM 

174 model (Girard et al. 2014 ) has the form 

\nn = £,+Bs, ( 13 ) 

175 where c, defines the terrain-following vertical coordinate, n denotes the hydrostatic pressure, 

1 76 B is a metric term prescribing the rate of flattening of the vertical coordinate with elevation, 

177 and 5 = ln(7T Jp re f) with n s being the hydrostatic pressure at the surface and p re f — 10 3 hPa 

1 78 is a reference pressure. The definition of this vertical coordinate follows the concept of the 

179 generalized hydrostatic-pressure-type hybrid coordinate proposed by Laprise ( 1992 ). Further 


9 



180 details regarding the log-hydrostatic-pressure-type vertical coordinate are provided by Girard 

181 et al. (2014) and Husain and Girard (2017). 

182 

183 The present study proposes to develop a GEM dynamical core where the vertical coordinate, 

184 given by Eq. (13), in the existing dynamical core is replaced by a height-based TFC. The tradi- 

185 tional formulation for height-based TFC can be expressed as 

C(z)=H (14) 
ZT Zs 

186 where z is the the true geometric height, zs and zt are the surface and the model top level heights, 
is? and H is a scaling constant. A more general form of Eq. (14) can be devised as 

z=A+Bzs, (15) 

i88 where A = (zt/H)£ and B — ( 1 — £ /H). Assigning H — zt implies zt = Cr and, as a result, Eq. 
iso (15) becomes 

z=£+Bzs, (16) 

190 which is similar to Eq. (13) in form. The vertical coordinate for the proposed dynamical core in 

191 the present study, however, is further generalized as 

z — C +BiZSL+B 2 (zs — Zsl), (17) 

192 which follows the concept of SLEVE (Smooth LEvel VErtical)-like coordinate system proposed 

193 by Schar et al. (2002), where zsl denotes the large-scale components of the orography. The vertical 

194 coordinate defined by Eq. (17) permits separate rates of flattening for the large and small-scale 

1 95 contributions of the orography on the terrain-following vertical coordinate with changing elevation 

196 through the metric terms B\ and Bi that are defined as 

b »=(MP <i8) 
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197 where r n — [hi,max ( r n,max r n,min)^k] ‘Mid A^ — (Cl Cfc)/(C 1 Cs)- The Values of F n ,min 3nd 

'as r n ma x together determine the rate of flattening of the vertical levels with increasing height. The 

199 subscript k of A indicates the model vertical level number. Furthermore, the value of k decreases 

200 with increasing height above the surface such that k= 1 indicates the top-most model level. 

201 

202 Henceforth, in this paper, the two dynamical cores with vertical coordinates based on log- 

203 hydrostatic-pressure and height are referred to as GEM-P and GEM-H, respectively. Different 

204 aspects of the GEM-P dynamical core, including the model formulation, discretization and 

205 numerical solution of the discretized problem along with the various modifications to the 

206 formulation over the past years, have been discussed in detail in the existing literature (Yeh 

207 et al. 2002 ; Qaddouri and Lee 2011 ; Girard et al. 2014 ; Husain and Girard 2017 ). The follow- 

208 ing subsections therefore only present the relevant details of the proposed GEM-H dynamical core. 

209 


210 c. GEM-H formulation 

211 The development of the GEM-H formulation requires further modifications to the system of 

212 equations (l)-( 5 ). First, virtual temperature, given by Eq. ( 12 ), is introduced in the system of 

213 equations along with an isothermal basic state temperature 7 ) such that T v = T' v + 7 ), where T', 

214 is the temperature deviation. The corresponding basic state pressure p* is defined hydrostatically 

215 as < 9 (lnp*) = —gdz/(RdE)- The equation of state, given by Eq. ( 11 ), is then used to eliminate 

216 density p as a prognostic variable followed by a transformation of the resulting equations from 

217 the geometric height coordinate to the terrain-following ^-coordinate. The vertical coordinate 

218 transformation leads to the replacements of the independent variables (jc,y,z) associated with the 

219 z-coordinate by (X, Y, 0 that are defined in the ^-coordinate. As a result, the system of equations 
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(l)-(5) is modified as follows: 
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dt 
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c pel Ty 


, Nl fdAnTy 
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(19) 

( 20 ) 

( 21 ) 

( 22 ) 

(23) 


phys 


where q = R c iT*\n(p //>*) is the pressure deviation from the basic-state pressure />*, t — yfr is the 
vertical motion with respect to the transformed ^-coordinate, — g 2 / ( c pd T *) is the square of the 
basic-state Brunt-Vaisala frequency and c 2 = R ( ]T X / (1 — K d ) is the square of the speed of sound. 
In Eq. (22), K is replaced by K d — R c i/c pc i as an approximation. It is also important to note that, 
the physical forcings associated with the modified continuity equation, given by Eq. (22), now 
includes the same diabatic heating term that appears in the thermodynamic equation, given by Eq. 
(23). 


232 In the above equations, the terms Jx, Jy and appears due to the vertical coordinate transfor- 

233 mation where Jx = j^, Jy — Jy and = Ji. It is however important to note that the coordinate 

234 transformation used to derive the Eqs. (19)—(23) is incomplete and only the following derivative 

235 operators have been transformed: 


236 


237 


d 

d 

Jx 

d 

dx 

~ dx 

h 

K 

d 

d 

Jy 

d 

dy 

~ dY 

h 
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d 

1 d 



dz 

hK 
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(24) 

(25) 

(26) 
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238 


d d d d-d 
Jt = Jt +U dX +V dY +Q dC' 


( 27 ) 


239 As a result, the original vertical velocity w has not been completely eliminated from the system. 
2« The system of equations in the ^-coordinate has its own vertical velocity in the form of t, = 

241 W’ whereas w = 4 ) remains in the system as a kinematic relation that needs to be dealt with 

242 explicitly. Charron et al. ( 2004 ) have demonstrated that such an approach is perfectly equivalent 

243 to a full coordinate transformation. As the treatment of advection in GEM is based on the semi- 

244 Lagrangian approach (Husain and Girard 2017 ), the kinematic relation defining w is also solved 

245 semi-Lagrangially. However, for convenience, the kinematic relation is modified as 

j t (z-Q + ^-w = 0, ( 28 ) 

246 where Eq. ( 28 ) along with the Eqs. ( 19 )—( 23 ) constitute the complete system of equations for the 

247 GEM-H formulation. 

248 One important aspect of any NWP model is how the effects of the physics forcings, as presented 

249 in the RHS of the Eqs. ( 19 )-( 23 ), are accounted for as the model equations are integrated at each 

25 0 time step. One way is to resolve the dynamical equations in the absence of any physics forcing 

25 1 and then modify the solution with the parameterized forcing as adjustments outside the dynamics 

252 step. Another possible approach is to compute the physics forcing and combine their effects with 

253 the nonlinear terms in a semi-implicit way within the dynamics step. This important aspect of 

254 the GEM model with particular focus on its impact pertaining to the GEM-H dynamical core is 

255 discussed in further details in section 4 . 

256 
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257 d. Spatial grid and discretization 

25s The objective of this project has been to implement the option of a dynamical core based on 

259 height-type vertical coordinate in addition to the existing pressure-type coordinate in the GEM 

260 model. The strategy has been to add the new coordinate option with minimal changes to the 

26 1 dynamical core and the rest of the GEM model source code. Therefore, the spatial grid struc- 

262 tures in GEM-H, both in the horizontal and the vertical, are kept the same as those in GEM-P, 

263 which implies a staggered Arakawa C grid (Arakawa 1988) in the horizontal and a staggered 

264 Chamey-Phillips grid (Charney and A.Phillips 1953) in the vertical. The horizontal and vertical 

265 grid structures are presented in Fig. 1. 

266 In addition to being similar to GEM-P for the limited-area model (LAM) grids, the global grid 

267 system is also kept unchanged in GEM-H, and is therefore, based on a Yin-Yang grid system 

268 (Qaddouri and Lee 2011). The Yin-Yang system combines two overlapping latitude-longitude 

269 LAM grids to form a global grid following the Schwarz method for non-matching domain 

270 decomposition (Qaddouri et al. 2008) and thus avoids pole-related singularity and convergence 

271 issues associated with a conventional global lat-lon grid. Further details on the Yin-Yang grid are 

272 provided by Qaddouri and Lee (2011). 

273 

274 e. Discretization in time 

275 The general form of an individual equation in the system comprised of Eqs. (19)—(23) and (28) 

276 can be expressed as 

dFj 

—77 ~ Gi = Pi (29) 

dt 

277 where F, denotes the advected quantity for an individual equation i within the system, G, is the 

278 associated dynamics source term with linear and nonlinear components, and P\ denotes the corre- 
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sponding physics forcing. Similarly to GEM-P, treating the advection terms in a semi-Lagrangian 
way and applying a two-time-level Crank-Nicholson temporal discretization leads to 


?d 


At 


l _±h G A 


1 -b, 


' G? = s c Pj 


(30) 


28 1 where At indicates the time-step length, and the superscripts A and D denote the arrival and depar- 

282 ture positions of the air parcels at the current time t and the previous time (/ — At), respectively. 

283 The integrals of the source terms G, for the different dynamical equations are approximated by 

284 trajectory averages. The parameter bj denotes the off-centering weight factor for the averaging 

285 of the dynamics source terms. When bj = 0, the averaging of the source term is fully centered, 

286 whereas bj > 0 implies additional weight placed on the implicit component of the source term. 

287 Historically, off-centering was implemented in GEM-P primarily to address spurious resonance 

288 originating from stationary orographic forcing (Rivest et al. 1994). However, it also suppresses 

289 computational noise and improves numerical stability. These other beneficial impacts have been 

290 found to be equally important in the current and previous implementations of GEM-P for the 

291 different operational NWP systems at ECCC. As the principal objective of this study is to have 

292 a GEM-H dynamical core that is equivalent to GEM-P for the different operational GEM-based 

293 NWP systems, off-centering has been retained in GEM-H. Also, following the latest implemen- 

294 tation of GEM-P (Husain and Girard 2017), a differential approach for off-centering has been 

295 adopted for GEM-H where bj varies depending on the dynamical equation denoted by the sub- 

296 script i. At present, the system of equations, given by Eqs. (19)-(23) and (28), are separated into 

297 three groups with an associated off-centering parameter for each group as follows: 


(i) b m for the horizontal momentum equations [Eqs. (19)—(20)], 


299 (ii) &/, for the continuity and thermodynamic equations [Eqs. (22)—(23)], and 


3 oo (iii) b n h for the vertical momentum and kinematic equations [Eqs. (21) and (28)]. 
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301 On the RHS of Eq. (30), the term P\ denotes the parameterized physics source term and the 

302 parameter s c indicates the mode of coupling between physics and dynamics. Depending on the 

303 chosen method for dynamics-physics coupling, the value of s c can be either 0 or 1. Also, the 

304 approach for dynamics-physics coupling determines how the contribution of P[ is accounted for in 

305 the model. Further discussions regarding the coupling of the parameterized physics forcing with 

306 the dynamical core is presented in section 4. 


so? f Trajectory calculations 

308 Semi-Lagrangian treatment of advection requires the solution of ki nematic displacement equa- 

309 tions of the form 


dX 

dt 


dY 

~ U '^t~ V ' 


d A = t 

dt ^ 


(31a) 

(31b) 


310 to determine the departure positions of the air parcels. In the context of GEM-P, Husain and Girard 

311 (2017) have shown that the consistency in the numerical discretizations between the dynamical and 

312 trajectory equations is fundamentally important for accurate solution of the advection problem. 

313 In order to be numerically consistent, similarly to the treatment of the dynamics source term in 

314 Eq. (29), the averaging of the velocities in Eq. (31) needs to be done using the trapezoidal 

315 rule. Furthermore, the interpolation scheme employed to determine the wind field at the departure 

3 1 6 positions for the trajectory calculations need to be the same as the one applied to determine the 
si? source terms in the dynamical equations at the departure positions. In the case of GEM-P, cubic 
sis interpolation is used for both the wind field and the dynamical source terms to achieve numerical 

3 1 9 consistency. Following the conclusions of Husain and Girard (2017), similar consistent trajectory 

320 calculation approach is adopted in GEM-H, i.e., trapezoidal rule for evaluating the integral of the 
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321 source term in Eq. (31) along with cubic interpolation to determine the wind field at the departure 

322 positions. 


323 g. The elliptic problem 


324 In order to solve the system of equations associated with the GEM-H formulation, each equation 

325 of the form (30) is rearranged to separate the linear and nonlinear components of the implicit part 

326 and is expressed as 

Li = Ri-Ni, (32) 

322 where ft = (Ff/x, - Gf) linear , N, = Ff'/x - G$ -L u and R, = FP/Xi + ftG? + s t ft with ft = 
328 ( 1 — ft)/(l + ft) and x ,• = Aftl + ft)/2. Eqs. (19)—(23) and (28) are rearranged as in Eq. (32), 


giving 




L u =-b 8 x q - Si JxJr 1 8 c q , 


n 


V 1 ^ b 

L v — - b 8yq — Sj/y-ft 8^q , 


Lw — ——b (Sj Jr 1 +Sd)5^ — g-ft 


^nh 


c 


L c = — (% -bln/H -b 8 x u-\ —— ftftcosft,) -bcft£ -ew^, 

% \G / COS0 


Lr = b(i * 


\ 

x h \T v c^/Tj 

Z-C 


L z = 


^?ih 


+ C -W, 


(33) 

(34) 

(35) 

(36) 

(37) 

(38) 


335 where the symbol ft denotes the finite difference operator along the /-direction, and the overline 

336 operator () J implies spatial averaging in the j-coordinate. For convenience of notation the terms 

337 4 an d -4 have been replaced by £ and /i, respectively. The corresponding nonlinear components 

C* o 

338 Nj associated with the discretized forms of the Eqs. (19)—(23) and (28) as well as the associated 

339 RHS terms R, are provided in Appendix A. The parameters Si and in the above equations denote 
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340 the choice of the solver for the dynamical core, where the subscripts i and d stand for iterative and 

341 direct respectively. Based on the selected solution approach, these parameters can be either 1 or 

342 0, and are mutually exclusive such that Sj = 1 — s^. Further discussion on the solution approaches 

343 is presented in Section 3. As shown by Cote and Staniforth (1988), the solution of the system 

344 of equations of the type (32) requires nonlinear iterations for convergence, where the nonlinear 

345 terms Nj are re-evaluated during each iteration using the latest values of the prognostic variables. 

346 Furthermore, Crank-Nicholson iterations are required, where the Rj terms are re-evaluated during 

347 each iteration at the departure positions calculated using the latest velocity estimates. At present, 

348 the GEM-based operational NWP systems at ECCC utilize two Crank-Nicholson iterations and 

349 within each Crank-Nicholson step two nonlinear iterations are carried out. As a result, irrespective 

35 0 of the solution approach, the solver is called four times during each dynamical time step. 

351 The discretized equations with the left hand sides (LHSs) given by Eqs. (33)—(38) are then 

352 reduced into a single elliptic boundary value (EBV) problem through elimination of variables, 

353 where the LHS of the final elliptic problem has the form 


L=8 X A + 


-8 Y [cos(j)B} + 8^C-eC- 


( —n— 

354 where A = I 8xq — Sj JxJ^ fyq J , B — I 8 Y q — $\J Y J^ 8^q J and C = F (s \J^ + Sd) 8^q — 

355 jiqt . Also, in Eq. (39), y = l/(c*T/,T,„) and F = 1 /(t m /z n h + N}zi,z m ). It is important to note 
sse that, with Si = 1, the terms A, B and C (with /i = 0) are simply the components of the gradients 

357 in the terrain-following coordinate system. The sequence of steps involved in deriving the EBV 

358 problem, i.e., the final form of L ", is provided in Appendix B. 
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359 


h. Initial and boundary conditions 


360 


361 


362 


363 


364 


365 


366 


As with the case of any initial value problem, in order to initiate integration in time, the GEM-H 
dynamical core requires initial values of all the prognostic variables. At present, ECCC’s opera¬ 
tional data assimilation system provides analyzed initial values for the horizontal wind components 
u and v, virtual temperature T v and surface pressure p s . The remaining prognostic variables w, £ 
and q are computed in a diagnostic manner at time t — 0 . The initial value of q is obtained from the 
analyzed surface pressure p s with a hydrostatic approximation. Substituting ^ = 0 in Eq. ( 21 ) 
gives 


1 dq _ T v ' 

T^dl =g Y v 


( 40 ) 


36 ? as a hydrostatic approximation. The value of q at the different model levels are then obtained by in- 

368 tegrating Eq. ( 40 ) where at the surface, due to the hydrostatic approximation, q s = R^T,. In (p s //?*). 

369 The initial value of £ is computed by assuming ^ = 0 at time t = 0 in the continuity equation. In 

370 the ^-coordinate, this takes the form 

d ( dz\ 1 d ( dz\ d ( f.dz\ 

dx{ pu 3s) + ^M ( C0S ^ v 5cJ + 5? H) = 1 °’ <41) 

371 which is then discretized to compute the initial value of £. Once t, is known, the initial value of w 

372 is obtained from its definition in the ^-coordinate, i.e., w=^ = uJx+vJy + C- 

373 The boundary conditions for the upper and lower boundaries are given by Cr = Cs — 0 . This 

374 implies that the vertical motion in the ^-coordinate vanishes at the surface and the model top which 

375 is flat. For LAM problems, GEM-H also requires lateral boundary conditions which are obtained 

376 from the driving fields. As the global Yin-Yang system is based on two interacting geometrically 
s?? identical LAM domains, it therefore similarly requires lateral boundary conditions. In this case, 

378 the boundary conditions for one sub-domain (Yin or Yang) depend on the solution in the other. 

379 Thus the solution of the global problem is obtained by iteratively solving the two sub-problems 
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380 separately and updating the values in the overlapping region until a certain convergence criteria is 

381 satisfied (Qaddouri and Lee 2011). 

382 3. Solution of the EBV problem 

383 The EBV problem to be resolved by GEM-H at each model time-step can be expressed as 

V^q + Mq — R, (42) 

384 by replacing l!" in Eq. (39) with (Vj. + M)q, where V| = (Sxx + ^ 0 <>Y (cos </>Sy)) is a discretized 

385 two-dimensional horizontal operator in the ^-coordinate, M contains all the remaining terms of 

386 L," that include the discretized difference and averaging operators in the horizontal and vertical 
38? dimensions, q is the unknown and R includes the explicit RHS terms as well as the implicit non- 

388 linear terms. It is important to note that the nonlinear terms in R require iterations for convergence, 

389 irrespective of the solution approach. At present, the GEM dynamical core uses two iterations for 

390 the sufficient convergence of the nonlinear terms and two iterations for trajectories. As a result, 

391 the solver is called into action four times during each dynamical step. 

392 Once Eq. (42) is solved to obtain the unknown q, the other prognostic variables are obtained 

393 through back substitution as presented in Appendix C. As has been mentioned earlier, two general 

394 approaches are available for solving the elliptic problem - direct and iterative. The selection of 

395 these approaches depends on which terms are included in M, and is determined by the values of 

396 the terms Si and s<i in Eqs. (33)-(37). 

397 a. The direct solver 

398 The direct solver works by first decoupling Eq. (42) in the vertical. It is achieved through the 

399 expansion of the unknown q and the RHS R in terms of the eigenvectors that diagonalizes the 

4 00 operator M (Qaddouri and Lee 2011). This is only possible when the operator M does not in- 
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401 elude contributions from the metric terms arising from the vertical coordinate transformation that 

402 involves horizontally variable coefficients imparting horizontal coupling. Therefore, the imple- 

403 mentation of direct solver in GEM-H requires a ‘simplified approach’ where all the metric terms 

404 of the relevant discretized equations are treated as nonlinear terms. This is achieved by setting 

405 S{i — 1. For N k number of vertical levels used in the model, vertical separation reduces Eq. (42) to 

406 a set of Nk independent horizontal Helmholtz problems of the form 

V 2 ^q + mq = R, (43) 

407 where q and R are the vertical projections of q and R, respectively, and m is the eigenvalue of the 

408 operator M. The horizontal solution of the algorithm then proceeds by expanding q and R in terms 

409 of the eigenvectors that diagonalize the X-component of the two-dimensional operator V|. For N, 
4 !o number of grid points along the longitude X, this leads to Nj independent tridiagonal problems of 
4 n Nj dimension for each model vertical level, where Nj denotes the number of grid points along the 

4 1 2 latitude Y. The total number of tridiagonal problems to solve is therefore Nk x Nj. Solution to the 

413 tridiagonal problems are computed by Gaussian elimination without pivoting, and afterwards, the 

4 1 4 final three-dimensional solution q is reconstituted (Qaddouri and Lee 2010). 

4 1 5 The ‘simplified approach’ has been primarily implemented in GEM-H to take advantage of the 

4 1 6 computational performance of the direct solver. The direct solver uses Fast Fourier Transform 
41? to compute the horizontal solution q and outperforms any iterative approach implemented at the 

4 1 8 ECCC by a substantial margin for the configurations of the various GEM-based NWP systems 

419 running operationally. The ‘simplified approach’ thus makes the application of GEM-H for such 

420 configurations feasible and keeps the option of a future replacement of the GEM-P core with 

421 GEM-H. The simplified approach, however, works as long as the vertical-horizontal coupling, im- 

422 parted through the metric terms, is not too significant so that these terms can be treated efficiently 
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423 through nonlinear iterations. This approach works unless the maximum terrain slope is not suffi- 

424 ciently steep (less than 30°) which is generally the case for most of the operational GEM-based 

425 NWP systems. However, with increasing spatial resolution, the slopes in grid-scale orography also 

426 increase, particularly over complex terrain, which leads to increased vertical-horizontal coupling 

427 and, at one point, makes the ‘simplified approach’ inapplicable. 

428 b. The iterative solver 

429 When all the metric terms in A and B (see Eq. 39) are included in M, the resulting vertical- 

430 horizontal coupling makes the problem non-separable. This is done by setting Sj = 1. In such a 

431 scenario, the three-dimensional equation of the form (42) can only be solved at each time step by 

432 using an iterative solver. The current iterative solver for the EBV problem in GEM-H is based 

433 on the flexible generalized minimal residual (FGMRES) method (Saad 1993; Qaddouri and Lee 

434 2010). 

435 The fully discretized system of equations of the form (42) can be further generalized as 

Ag = R (44) 

436 where coefficient matrix A contains the discretized operator (V| +M). The FGMRES method 

437 approximates the solution in a Krylov sub-space of small dimension by minimizing the Euclidean 

438 norm of its residual. A major advantage of such a method is that instead of explicitly generating 

439 the coefficient matrix A, one only needs to compute the vector resulting from the action of the 

440 underlying operator (V| + M) on the vector q. Efficient functioning of such an iterative solver, 

441 however, requires a pre-conditioner. At present, the pre-conditioner is based on the block Jacobi 

442 iteration for the EBV in Eq. (42), where all the metric terms in M are absent. This pre-conditioner 

443 improves the convergence rate of the FGMRES solver. However, the time of execution is still high 
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444 compared to the fast direct solver. Significant research is currently underway at ECCC to devise 

445 iterative solvers that are competitive with the direct approach and will work for both GEM-P and 

446 GEM-H. At present, the current implementation of the iterative solver in GEM-H—although not 

447 as efficient—provides the necessary reference for the direct solver approach. 

448 4. Dynamics-physics coupling 

449 Along with the dynamical core, parameterization of the subgrid-scale physical processes consti- 

450 tutes the other fundamental component of any NWP model. Coupling between the dynamical core 

45 1 and the parameterized subgrid-scale physical processes is of critical importance. How to devise 

452 the most appropriate coupling strategy is still an unsettled question (Beljaars et al. 2018). This 

453 issue is being actively studied at ECCC. However, it is not the objective of this study to delve 

454 deep into the fundamental questions around dynamics-physics coupling. Rather, in this section, 

455 the issue of coupling the RPN Physics Package with the GEM dynamical core is discussed in order 

456 to determine which approach is the most feasible for GEM-H among the options that are available 

457 for GEM-P. 

458 It is important to note that, during every model time step in GEM, the RPN Physics Package 

459 is executed after the dynamical equations have been resolved by the dynamical core and thus the 

460 physics schemes utilize the outputs of the dynamics step as inputs. However, irrespective of the 

461 vertical coordinate used in the dynamical core, the physics schemes use a traditional c-coordinate 

462 in the vertical, defined as 

K 

<x=- (45) 

K s 

463 where n is the hydrostatic pressure. Also, the physics schemes work within a one-dimensional 

464 configuration where each processor core only has access to the vertical structure of the meteoro- 

465 logical fields associated with a single horizontal grid point. The various physical processes are 
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466 parameterized sequentially where the tendencies estimated by one parameterization scheme af- 

467 fects the ones that follow. The parameterization sequence during each physics step initiates with 

468 the radiation scheme and is followed by the parameterizations of the surface processes, gravity 

469 wave drag, boundary layer turbulence, convection and grid-scale condensation. At the end of the 
4/0 physics step, the tendencies from the different physical parameterization schemes are aggregated 

471 to compute the grid-scale tendencies for the wind components, temperature, water vapor and the 

472 hydrometeors. 


473 a. Different coupling approaches in GEM 

474 Particularly in the context of GEM-P, two approaches are presently available to couple the RPN 

475 Physics Package with the GEM dynamical core. A brief discussion on these methods will be 

476 helpful in establishing the rationale behind the approach selected for the GEM-H dynamical core. 


(i) Split method: As has been mentioned earlier, s c determines the mode of coupling between 
dynamics and physics. If s c = 0 then the dynamical equations are resolved in the absence of 
any physical forcing and at the end of the dynamics step their contributions are incorporated 
as adjustments in the so called ‘split mode’. In the absence of physics forcing, Eq. (30) 
becomes 


pA* _ pD 
l _ j 

At 


1 + bi r A* 
~ a ‘ 


i-b, D 

— tr Gi 


0 


(46) 


where F/ 4 * is the interim solution of the dynamics step. Once Eq. (46) is resolved, the physics 
source term is then applied as grid point adjustments as follows 


(SF) plr ,, = F^-Fp = AtP,. 


(47) 
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Thus, in the split method, the dynamics step predicts an inviscid and adiabatic solution that 
is modified through adjustments attributable to the parameterized physics forcings in order 
to obtain the complete solution at the end of each model time step. 


487 


488 


489 


490 


(ii) Explicit method: The second option for dynamics-physics coupling treats the physics source 
terms explicitly by setting s c = 0 and replacing / J ( by Pf. This method is referred to as the 
‘explicit method’ and moves by solving equations of the following form at each model time 
step 


At 


D 


1 + bi r A 

—G, 


1 -bi 


Gf = Pi 


i D 


(48) 


In this approach, physics tendency Pj from the previous time step is combined with the RHS 
term Rj, followed by the determination of Rf at the departure positions in a semi-Lagrangian 
way. In other words, the explicit method works by directly incorporating the impact of 
physics forcings as tendencies into the discretized dynamical equations. 


495 b. Coupling in GEM-H 

496 Although both of the aforementioned coupling methods are available for GEM-P, it is important 

497 to note that all the operational NWP systems based on GEM-P at present utilize the split method 

498 for dynamics-physics coupling. Nevertheless, there exists strong concern about the split method 

499 in general, and a brief discussion highlighting the pertinent issues will be helpful. 

soo In the context of model formulations for both GEM-P and GEM-H, the term Fj does not nec- 
5 oi essarily coincide with the prognostic model variables. For example, in the case of GEM-P in its 
soo hydrostatic mode, as presented by Girard et al. (2014): 

E = jw, v ,Bs + ln^l + lnp- K d Bs^j , (49a) 

P = {^( M ’ v ’ ln P’ lnrv )} » (49b) 
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503 whereas the prognostic variables are u, v, s, £, and T v . Therefore, in the actual model implemen- 

504 tation of the split method, the prediction from the dynamics step is utilized to compute the interim 

505 state of the prognostic variables and adjustments are then applied to these variables at the end of 


506 the physics step. It is important to note that in this case, the only adjustments that have been found 
so? to not result in any issue of major concern are ( 8u) p hy S , (Sv) p f iys and ( 8T v ) p h ys . Also in GEM-P, 
508 adjustments are required for density, as in effect 


( d In p \ 

V dt ) phys 


d_ 

dt 


ln(l + r w ). 


(50) 


509 Although water vapor and hydrometeors are updated through physical parameterizations, the only 


5 10 variable that could be updated in the split mode appears to be s through a surface pressure adjust- 

511 ment 8k s that may be computed as 


8n s — 8 ln(l + r w )dn, 


5 i 2 which takes into account the net inflow/outflow of mass through the Earth’s surface at every model 
sis grid point. Unfortunately, the vertical distribution of this change in mass through water vapor and 

514 precipitation fluxes cannot be correctly accounted for in the split mode, and is found to produce 

5 1 5 considerable noise in the wind forecast. 

5 1 6 Furthermore, in the context of three-time-level discretization, Caya et al. (1998) have shown 
si? that the split method can lead to erroneous results for long time steps that are permissible with 
sis the semi-implicit semi-Lagrangian models. Similar conclusions were drawn from a theoretical 

5 1 9 analysis for two-time-level schemes by Staniforth et al. (2002). Figure 2a shows the geopotential 

520 height contours at 400 hPa from a 72-h global forecast with ECCC’s 25-km resolution Global De- 

521 terministic Prediction System (GDPS) using the GEM-P dynamical core. The results correspond 

522 to split method for dynamics-physics coupling. Although the distribution of geopotential height 

523 for the meso and large scales does not reveal any issue of immediate concern, when one looks at 
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524 the smallest scales, i.e., scales of about a few grid lengths, some spurious computational noise is 

525 visible (see Fig. 2 b). Historically, the operational NWP systems at ECCC have been utilizing spa- 

526 tial filters over the model-predicted meteorological fields of interest, like the geopotential height, 

527 to smooth out any computational noise in the model outputs. As a result, the kind of computa- 

528 tional noise shown in Fig. ( 2 b) has not been troublesome for the meteorologists using ECCC’s 

529 operational NWP outputs. Nevertheless, the computational noise associated with the split method 

530 remains as a concern. However, as shown in Fig. 2 c, when the split method for coupling is re- 

531 placed by the explicit method, the noise in the geopotential height disappears. It should be noted 

532 that, although explicit coupling can impose stability limitations in terms of the acceptable length 

533 of time steps, ECCC’s operational NWP system configurations are found to function with explicit 

534 coupling without requiring any adjustments to the time-step lengths. 

535 Similarly to GEM-P, the F t terms do not coincide with the prognostic variables for all of the 
sse GEM-H model equations, where 


r ) <7 . , , , T v Q 

F=^,v,w, 5+ 1„7 c ,1„--— r 


phys 


p = { J t \ u,v, W ,m P T v MT v 


( 52 a) 

( 52 b) 


537 Particularly, the presence of the physics forcing term ( dln f Tv ) in the modified continuity 

538 equation ( 22 ) poses an additional challenge for the GEM-H formulation, as far as the split method 

539 is concerned. If one attempts to apply this tendency as a hydrostatic adjustment to pressure in 

540 the split mode then of course its effect is not limited to the continuity equation alone. Rather, 

541 applying the adjustments to pressure, i.e. changing q in the case of GEM-H, also affects the 

542 thermodynamic equation. Such an adjustment leads to over compensation and is found to result 

543 in spurious bias in the temperature in upper troposphere and stratosphere, which is unacceptable 

544 (not shown). Also, the jet-level wind is found to be adversely affected. As a result, in GEM-H the 
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545 physics contributions are accommodated through the explicit method by including them as explicit 

546 grid-scale tendencies in the RHS of the discretized dynamical equations. For a fair comparison 

547 between the new and the existing dynamical core, the results for both the GEM-P and GEM-H 

548 cores presented in the rest of this paper are obtained with explicit dynamics-physics coupling. It is 

549 also important to mention that, even though parameterization of physical processes like boundary- 

550 layer turbulence can modify vertical motion w, at present the impact of physics on w is neglected. 

551 This is the case for both GEM-P and GEM-H. 

552 A hybrid ‘split-explicit method’ has also been tested with GEM-H, where the physics contri- 

553 butions to the thermodynamic and horizontal momentum equations are accounted for through the 

554 ‘split method’ while the contributions to the continuity equation is accommodated using the ‘ex- 

555 plicit method’. Such a hybrid approach with GEM-H produces results that are equivalent to those 

556 obtained with the ‘split method’ for GEM-P. However, questions remain about the numerical con- 

557 sistency of the hybrid split-explicit method. 

558 It is important to mention that the current implementation of the RPN Physics Package, when 

559 coupled with the GEM-P (GEM-H) dynamical core through the explicit method, leads to some 

560 deterioration in temperature bias in the upper troposphere compared to the split (split-explicit) 

56 1 method. This implies that further research is necessary to have a more consistent dynamics- 

562 physics coupling with the explicit method. Particularly, the grid-scale condensation scheme in 

563 GEM, for 10 km or coarser horizontal resolutions, has been found to exhibit large sensitivity with 

564 the explicit method which leads to under-prediction of clouds (R. McTaggart-Cowan, ECCC, per- 

565 sonal communication). This implies that some process-specific adjustments in the computation 
see of the relevant physics tendencies may be required to improve the overall dynamics-physics cou¬ 
se? pling. The challenges imposed by the process-specific issues are also being explored by other 
sea operational NWP centers (Beljaars et al. 2018 ). Currently, work is underway at ECCC to explore 
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569 the various issues within the coupling interface as well as the parameterizations of the different 

570 physical processes to improve the dynamics-physics coupling in general. Further discussion on 

571 this issue, however, is beyond the scope of this paper. 

572 5 . Evaluation of GEM-H 

573 One of the most important objectives of this study has been to develop a dynamical core for 

574 the GEM model that utilizes a height-based TFC and is capable of producing predictions that 

575 are equivalent to the results obtained by the pressure-based dynamical core. In order to evaluate 

576 the consistency and performance of the new dynamical core, a number of numerical experiments 

577 covering a wide range of scales—ranging from microscales to the meso and synoptic scales—have 

578 been carried out. These include two-dimensional theoretical test cases involving bubble convection 

579 (Robert 1993 ) and nonhydrostatic mountain waves (Schar et al. 2002 ) as well as three-dimensional 

580 global NWR The two-dimensional test cases selected in this paper have become ubiquitous tools in 
sal testing the consistency and performance of nonhydrostatic dynamical cores. Also, the availability 

582 of the GEM-P core provides the opportunity to have reference solutions for all these cases. 

583 a. Robert’s bubble convection case 

584 Robert ( 1993 ) presented a two-dimensional theoretical case involving the evolution of a warm 

585 bubble within a dry isentropic atmosphere. Initially, the bubble has a diameter of 500 m and is 

586 placed 10 m above a flat surface within a 1 km x 1 km computational domain, and has a uniform 

587 potential temperature of 30 . 5 °C. Also, the basic-state atmosphere is at rest under a hydrostatic 

588 equilibrium with an isentropic basic-state temperature of 30 °C. As the bubble has a potential 

589 temperature excess of 0 . 5 °C compared to the surrounding atmosphere, it rises due to the action of 

590 the buoyancy force. The absence of any orographic variation makes this experiment an excellent 
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591 benchmark to test the functioning of advection and buoyancy within a dynamical core during the 

592 early stages of its development. 

593 The numerical experiment for bubble convection is carried out with a spatial grid resolution of 

594 10 m and a time step of 5 s. No explicit numerical diffusion is used. The resulting evolution 

595 of the bubble, in terms of its potential temperature distribution at two different times (7 min and 

596 10 min), is presented in Fig. 3 for both GEM-P and GEM-H. The bubbles predicted by the two 

597 cores initially deform into a somewhat mushroom-like shape (see at r = 7 min) and then are de- 

598 formed further (at t — 10 min). Overall, the predictions from GEM-P and GEM-H are equivalent 

599 for the entire range of scales - from the large to the smallest scales. Such a good resemblance be- 

6 00 tween the two predictions imply negligible impact of the choice of the vertical coordinate and the 

6 01 other modifications in model formulations in the absence of any orographic variation at the model 

602 surface. It also indicates that the representation of the advection and buoyancy effects are compa- 

603 rable between the two GEM cores. It is important to note that, due to considerable differences in 

604 the model formulations and spatiotemporal discretizations, it is difficult to compare the evolution 
60 s of the bubbles between two completely separate models in a quantitative manner. Only qualita- 
606 tive comparisons are feasible. Therefore, the lesser resemblance between the results from Robert 
60 ? (1993) and the GEM dynamical cores are not unusual. Although the predictions from the two GEM 
60 s cores have some large-scale resemblance to the results presented by Robert (1993), significant dif- 

609 ferences appear at the upper half of the bubble - particularly at t — 10 min. However, the upper 

6 10 structure of the bubble compares better with the predictions by Smolarkiewicz and Pudykiewicz 
on (1992). Also, because of the implicit dissipation associated with the semi-Lagrangian approach, 

612 the GEM solution does not suffer from computational noise like models based on Eulerian ad- 

6 1 3 vection (Juang 2000). Overall, as GEM-P is being used operationally at ECCC, the resemblance 

6 1 4 between GEM-H and GEM-P is of more significant importance, as it confirms consistency of the 


30 



615 GEM-H formulation and a neutral impact of the vertical coordinate modification on buoyancy and 

616 advection. 

si? b. Schar’s mountain case 

618 Schar et al. (2002) presented a linear two-dimensional theoretical test case of mountain waves 

619 which is an excellent benchmark for verifying nonhydro static dynamical cores, particularly in 

620 determining the presence of possible inconsistencies in the numerical details (Husain and Girard 

621 2017; Melvin et al. 2010; Girard et al. 2005; Klemp et al. 2003). The bottom boundary profile of 

622 the idealized mountain for this case is defined by 

z s = zjoe~W a ' ) cos 2 (nx/l x ), (53) 

623 where zo=250 m, l x =4 km, a =5 km, and n is the conventional mathematical constant. The upstream 

624 flow conditions are given by uniform upstream wind U=10 m s , upstream surface temperature 

625 T surf = 288 K, upstream surface pressure po= 1000 hPa, and a constant Brunt-Vaisala frequency 

626 A*=0.0 1 s' 1 . All other conditions for the simulations of this test case is similar to those presented 

627 by Husain and Girard (2017). 

628 One major advantage of Schar’s mountain case is the availability of a steady-state analytical 

629 solution of the corresponding to the linearized problem as a reference (Schar et al. 2002). The 

630 simulated quasi-steady vertical velocity obtained after 4 hours of integration with both the GEM-P 

631 and GEM-H dynamical cores are presented in Fig. 4. The analytical solution of the problem, as 

632 presented by Schar et al. (2002), is a combination of rapidly decaying small-scale nonhydrostatic 

633 mountain waves close to the surface and large-scale hydrostatic waves extending to much higher 

634 altitudes. As can be seen in Fig. 4, the solutions with both types of vertical coordinates generate 

635 these two regimes of the mountain waves and are very similar. Results shown in this figure repre- 
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636 sent simulations that have been carried out without any off-centering in the discretized dynamical 

637 equations and with consistent trajectory calculations, i.e., integration of the wind field in the tra- 

638 jectory equations is based on the trapezoidal rule while the wind field at the departure positions 

639 is obtained with cubic interpolation. Inconsistent trajectory calculations were found to produce 
6 M similar distortion in the large-scale hydrostatic waves with GEM-H (not shown) as has been found 

641 for GEM-P earlier by Husain and Girard (2017). Although the GEM-H results presented here cor- 

642 responds to the direct solver based on the simplified approach, results with the iterative solver is 

643 found to be almost identical (not shown). Figure 5 reveals that off-centering in the discretized dy- 

644 namical equations leads to distortions in the vertical velocity distribution, irrespective of the type 

645 of the vertical coordinate. The solutions presented here are obtained with uniform off-centering 

646 involving b m = bj, — b n h = 0.2, which are the standard values used in the current GEM-based 

647 operational NWP systems. The results correspond to At=32 s, for which the maximum Courant 

648 number is approximately 0.76. Reducing the time step to even 4 s is unable to remove these dis- 

649 tortions. Also, reducing the level of off-centering is found to reduce the level of distortion in the 

650 mountain waves, but the distortions are only completely eliminated when b m = /?/, = b,,/, — 0 (not 

651 shown). This conforms to the conclusions drawn by Husain and Girard (2017) in the context of 

652 GEM-P. Furthermore, as has been shown by Husain and Girard (2017), consistent trajectory calcu- 

653 lations in the presence of off-centering necessitates off-centered averaging applied to the integrals 

654 of the source term on the RHS of Eq. (31). With uniform off-centering applied to the discretized 
ess dynamical equations, the discretized trajectory equations also require uniform off-centering of the 

656 same degree. Figure 6 reveals that in the presence of consistent off-centering in the trajectory 

657 calculations, the distortions in the vertical velocity distribution are eliminated for the both dynam- 

658 ical cores. In the presence of differential off-centering, i.e., with different values of b /, and /?„/, 

659 for hydrostatic and nonhydro static contributions in the system of dynamical equations, a similar 
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660 differential approach is required for off-centering in the discretized trajectory equations. As has 

661 been shown for GEM-P by Husain and Girard (2017), in order to achieve numerical consistency, 

662 the off-centering in the discretized source terms in Eq. (31a) and (31b) need to be equal to the 

663 values of b /, and b n i t used in the discretized dynamical equations, respectively. 

664 c. Global deterministic prediction 

ees A series of 5-day global forecasts, with 25 km horizontal grid spacing, has been carried out 
eee covering winter and summer periods for the Northern Hemisphere to compare the predictions 
667 by GEM-H and GEM-P from NWP standpoint. Each seasonal period includes 44 cases where 
ees the initialization between two consecutive cases are apart by 36 hours. The first summer and 
669 winter cases start at 0000 UTC of 25 June 2014 and 19 December 2014, respectively. The global 
e/o predictions have been obtained with uniform off-centering in the discretized dynamical equations 

671 ( b m = bj, = b n h = 0.2). Husain and Girard (2017) have shown that inconsistent off-centering has 

672 negligible effects for three-dimensional NWP applications. Therefore, the global forecasts are 

673 carried out without any off-centering in the discretized trajectory equations. Furthermore, as has 

674 been mentioned in section 4, physics is coupled with dynamics through the explicit method for 

675 both GEM-H and GEM-P for all the tests presented in this section. 

676 First, comparisons are made in the spectral space by comparing the variance spectra associated 
e?? with the different meteorologically important fields. In order to compute the spectral variance of 

678 any meteorological field, it is first interpolated from the Yin-Yang grid to a global Gaussian grid. 

679 Afterwards, variance spectra of the field are calculated by decomposing the field at a given pressure 

680 level using the spherical harmonics. Figure 7 shows the spectral variance of the geopotential height 

681 and temperature fields for 120-h forecasts at three different pressure levels for an average over 10 

682 cases for the winter period. The results show spectral similarity between GEM-P and GEM-H 
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683 for the entire range of scales resolved by the global model - from synoptic to mesoscales. The 

684 spectra of kinetic energy and vertical velocity are presented in Fig. 8. The spectral slope of kinetic 
ess energy is critical for accurate representation of atmospheric dynamics, and as can be seen in this 
ese figure, both GEM-H and GEM-P have the same spectral slope at the synoptic and mesoscales. The 
68 ? vertical motion is also important, particularly for physical processes like convection. Fig. 8 shows 
ess close spectral similarity between the vertical motions from the two dynamical cores. This implies 

689 that changing the vertical coordinate has negligible sensitivity to the extremely important physical 

690 process like convection and convection-driven precipitation. Also, the comparisons in the spectral 

691 space confirms that the height-based TFC does not lead to any spurious noise or damping in the 

692 meteorological fields for any model-resolved length scale. 

693 Objective forecast scores are computed by comparing the model predictions against radiosonde 

694 observations at different pressure levels. The evaluation is based on the bias and standard deviation 

695 of error (SDE) in model predictions for the individual cases as well as for the average of the 44 

696 cases covering each seasonal period. Figure 9 presents the vertical profile of error in the predic- 
69? tions from GEM-P (blue) and GEM-H (red) for the winter period. These figures represent global 

698 average scores of 120-hour forecast from 44 winter cases for zonal wind (UU), wind speed (UV), 

699 geopotential height (GZ), and temperature (TT). An important thing to note while reading this 

700 figure is the presence of the statistical confidence scores at the different pressure levels along the 

701 left and right vertical axes of the individual subplots for bias and SDE, respectively. A confidence 

702 value (in %) shaded in blue (red) color implies statistically significant improvement obtained with 

703 the GEM-P (GEM-H) core with respect to the other. The confidence score for the average of SDE 

704 and bias are estimated by applying the F-test and t- test, respectively. Figure 9 reveals that although 

705 there are small differences in the bias for the average of the 44 winter cases, there is no statistically 

706 significant difference in the SDE. When tested in the absence of physics forcings, no statistically 


34 



707 significant difference is found between the two dynamical cores in either bias or SDE (not shown). 
70 s The meteorological fields are interpolated from the TFC of the dynamics to the c-coordinate for 

709 physics through vertical interpolation which can lead to small differences in the vertical for the 

710 different definitions of the TFC. Physical parameterizations can be sensitive to the position of the 
7 n vertical levels and is apparently responsible for the small bias differences shown in Fig. 9. Even 

712 though small differences in bias are present, the objective scores from GEM-P and GEM-H can be 

713 safely assumed to be equivalent as a whole. During the summer period, the two dynamical cores 

714 have also been found to be similarly equivalent (not shown). 

715 The geopotential height at 1000 hPa in Fig. 9 shows a negative bias for both dynamical cores. 

7 1 6 This indicates a loss of mass conservation, which is a consequence of the non-conservative nature 

717 of semi-Lagrangian advection. This can be improved by introducing a simple global mass fixer 

7 1 8 that works by conserving the global mean surface pressure after each dynamics step in the model. 

719 The scores in the presence of a global mass fixer is shown in Fig. 10, where the bias at the lowest 

720 model level is eliminated for both dynamical cores. The overall scores for the two cores are again, 

721 as expected, found to be equivalent in the presence of a global mass fixer. 

722 Overall, the results presented in the Figures 3-10 clearly demonstrate that the implementation 

723 of GEM-H as presented in this paper produces results that are equivalent to the existing GEM-P 

724 dynamical core. 

725 6. Summary 

726 A newly-developed dynamical core for ECCC’s GEM model with a height-based terrain- 

727 following vertical coordinate has been presented. With increasing focus on three-dimensional 

728 iterative solvers at ECCC driven by the limitations of the operational direct solver as well as the 

729 strong numerical instability induced by steep-orography for sub-kilometer resolution NWP, a dy- 
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730 namical core with height-based TCF is expected to be better placed to address the future NWP 

731 challenges at ECCC. The principal objective of this paper is to provide information pertaining to 

732 the different aspects of the new height-based dynamical core including changes to the model for- 

733 mulation, discretizations, solvers for the discretized problem and the strategy for coupling the new 

734 core with the RPN physics package. Another important objective is to demonstrate that the new 

735 GEM-H core is capable of making meteorological predictions that are equivalent to the existing 

736 GEM-P dynamical core, which is based on a log-hydrostatic-pressure-type vertical coordinate. 

737 Numerical experiments have been conducted throughout the different stages of GEM-H devel- 

738 opment. Initially, the bubble convection test revealed that the advection and the buoyancy effects in 

739 GEM-H are accurately represented and are producing results that are equivalent to GEM-P. When 

740 tested for the idealized Schar’s mountain case, the nonhydro static and hydrostatic components of 

741 the mountain waves predicted by GEM-H are found to be very close to the GEM-P predictions 

742 as well as the analytical solution. The dynamics source terms in GEM are averaged over the 

743 air parcel trajectories using the trapezoidal method and the calculation of the RHS terms in the 

744 dynamical equations are carried out using cubic interpolation. Although it has not been shown 

745 explicitly, similarly to GEM-P, in the absence of any off-centering in the discretized dynamical 

746 equations, numerical consistency in GEM-H requires a trapezoidal averaging of the source terms 

747 in the discretized trajectory equations along with a cubic interpolation for the wind fields at the de- 

748 parture positions. Furthermore, in the presence of off-centering, the Schar’s mountain case shows 

749 that the discretized sources terms in the trajectory equations also necessitate off-centered averag- 

75 0 ing for the sake of consistent numerics in both GEM-H and GEM-P. In general, the knowledge 

75 1 acquired over years regarding the different numerical aspects of the GEM-P dynamical core is 

752 proven to be equally applicable to the case of GEM-H. Comparisons between GEM-H and GEM- 

753 P for global deterministic predictions are also presented. The results are found to be equivalent in 
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754 the spectral space confirming that GEM-H does not produce any spurious noise or damping over 

755 the model-resolved scales. When compared against upper-air radiosonde observations, except for 

756 small differences in bias, GEM-H and GEM-P are found to produce equivalent results. 

757 The rationale behind the choice of the solution approach for the discretized EBV resulting from 

758 the GEM-H formulation as well as the method for dynamics-physics coupling is discussed in 

759 detail. Although the general structure of the EBV originating form the GEM-H discretized sys- 

760 tern of equations require an iterative approach, a simplified approach has been devised where 

761 the horizontally-variable metric terms—attributable to the vertical-coordinate transformation—are 

762 coupled with the nonlinear terms and are treated with nonlinear iterations. This makes the EBV 

763 vertically separable and allows the use of a direct solver which is computationally very efficient 

764 for the currently operational NWP system configurations. A three-dimensional iterative solver 

765 based on FGMRES is also developed for situations when the simplified approach is not feasible. 

766 The fact that the GEM-H core can utilize the direct solver approach for global and regional scale 

767 model resolutions, eliminates the concern of computational efficiency as far as its implementation 

768 in the current and near-future plans for operational NWP systems at ECCC is concerned. 

769 Improving any NWP model is a continuous process. As a stable GEM-H dynamical core is 

770 now developed, a number of other relevant issues are currently being studied. The objective 

771 is to improve the GEM model in general and the GEM-H dynamical core in particular. One 

772 of the most important short-term goal in this regard is to devise a more numerically consistent 

773 and accurate coupling between RPN Physics and GEM dynamics which will benefit both the 

774 dynamical cores. Also, extensive research is being conducted to develop highly optimized 

775 three-dimensional iterative solvers that can be competitive against the operationally-used direct 

776 solver while scaling better for very large number of processor cores for the future generations 

777 of massively parallel supercomputers. Currently, the Yin-Yang system uses the Schwarz method 
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778 where the global solution is produced by iteratively solving two elliptic sub-problems for the 

779 two sub-domains (Yin and Yang) separately and updating the solutions in the overlapping 

780 regions until a certain convergence criteria is satisfied. One promising solution to reduce the 

781 execution time of the iterative solver for the Yin-Yang system is to remove the Schwarz iterations 

782 and to solve the two elliptic sub-problems as one by using FGMRES (Zerroukat and Allen 

783 2012). Also, pre-conditioners based on other types of methods, e.g., incomplete factorization, 

784 block Gauss-Siedel or multigrid method, could be used in order to improve the convergence 

785 rate of the FGMRES solver. Finally, on the GEM-H front, another important short-term 

786 objective is to improve its numerical stability over steep orography by implementing more ac- 

787 curate numerical approximation of the horizontal gradients in the discretized dynamical equations. 


789 Acknowledgments. The authors would like to sincerely thank the members of the Numerical 

790 Methods Research Group at RPN for all their thoughtful comments and suggestions during the 

791 course of GEM-H development. The authors would like to particularly thank Michel Desgagne, 

792 Stephane Gaudreault, Rabah Aider and Vivian Lee for their contributions during the implemen- 

793 tation of the GEM-H source code. Also, comments from the internal reviewer, Dr. Christopher 

794 Subich, have helped to considerably improve the overall presentation of the paper. 


APPENDIX A 


Nonlinear and RHS components of the discretized equations 


The nonlinear components of (19)—(23) and (28), associated with linear components (33)—(38), 
are given by 


Nu = -\ fAy + - 1W -cy -s, 


(Al) 
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The corresponding RHS terms, in the absence physics forcing, take the following forms 
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8,o APPENDIX B 

8 „ Deriving the EBV problem 

8 i 2 The LHS of the EBV problem to be solved at every time step in GEM-H is derived by manipulating 

8.3 the discretized system of equations of the form (32). The sequence of operations to derive the EBV 

8.4 are provided below in terms of its linear components: 

, 1 ln/r\ 
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Similar manipulations are also applied to the nonlinear and RHS components of the discretized 
equations to obtain R c , n' c , R z , N z , R w , N w , R t , n' t , R t , Nj, R." c , n", r”. and n" ( '. 


APPENDIX C 


Back-substitution to obtain the other prognostic variables 

The solution of the EBV computes the unkown q. The rest of the variables are then calculated 
from the discretized dynamical equations through back-substitution in the following sequence: 
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931 Fig. 10 . Same as in Fig. 9, but with a global mass fixer in the simulations for both GEM-P and 
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Fig. 1. a) Vertical Chamey-Phillips grid, b) Horizontal Arakawa C grid. 
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933 Fig. 2. a) Geopotential height (dam) at 400 hPa after 72 hours for a 25-km global forecast initiated at 1200 

934 UTC of 25 January 2015 obtained with GEM-P using the ‘split method’ for dynamics-physics coupling (contour 

935 intervals of 6 dam), b) The enlarged view over a small section of the global domain (contour intervals of 0.5 

936 dam), c) Same as in Fig. 2b, but with ‘explicit method’ for the dynamics-physics coupling. 
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a) 


GEM-P (f = 7 min) 


b) 


GEM-H (t= 7 min) 



937 Fig. 3. Potential temperature (°C) distribution for the rising bubble with: a) GEM-P at t — 7 min, b) GEM-FI 

938 at t = 7 min, c) GEM-P at t — 10 min, and d) GEM-FI at t — 10 min. 
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941 Fig. 5 . Same as in Fig. 4 , but with uniform off-centering (b m — b /, = = 0 . 2 ) used in the discretized 

942 dynamical equations. The absence of off-centering in the discretized trajectory equations leads to inconsistent 

943 trajectory calculations and distortions in the mountain waves. 
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a) GEM-P 



Fig. 6 . Same as in Fig. 4 , but with consistent off-center 
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Fig. 7. Spectral variance of (left column) geopotential height (GZ) and (right column) temperature (TT) for 


M 7 120-h global forecasts with 25 km horizontal grid spacing. Results are presented for three different pressure 


948 levels (top: 100 hPa, middle: 500 hPa, bottom: 850 hPa) that are obtained by averaging the spectra for 10 


949 Northern Hemisphere winter forecasts. 
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Fig. 8 . Spectral variance of (left column) kinetic energy (KE) and (right column) vertical velocity (WW) for 


951 120-h global forecasts. All other conditions as in Fig. 7. 
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952 Fig. 9. Comparison of 120-h 25-km GDPS forecasts obtained with GEM-P (blue) and GEM-F1 (red) dy- 

953 namical cores against radiosonde observations for zonal wind (UU), wind speed (UV), temperature (TT), and 

954 geopotential height (GZ). The dashed and solid lines, respectively, indicate bias and standard deviation of error 

955 (SDE). The scores are obtained by averaging over 44 Northern Flemisphere winter cases. The red and blue 

956 shaded numbers along the left (right) vertical axes within each panel indicate the confidence in percentage in the 

957 statistically significant improvements in bias (SDE) for the dynamical core associated with each color. Signifi- 

9 58 cance for bias and SDE are computed using t- and F-tcst, respectively. 
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959 Fig. 10. Same as in Fig. 9, but with a global mass fixer in the simulations for both GEM-P and GEM-FI to 
geo improve mass conservation. 
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